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ABSTRACT 


The  numerical  experiments,  carried  out  through  the  use  of  a  pressure-velocity 
coupled  method  to  solve  the  Favre  Averaged  Navier-Stokes  equations,  on  steady  and 
sinusoidally  oscillating  flows  at  five  different  Keluegan-Carpenter  numbers,  and 
three  periodicity  levels  are  described.  A  second-order  in  time,  second-order  in 
space,  second-level  predictor-corrector  finite  difference  scheme  has  been  used.  The 
solutions  were  solved  by  the  CFD-ACE  program  from  the  CFD  Research 
Corporation.  The  analysis  has  produced  in-line  force  coefficients  comparable  to  those 
obtained  experimentally  for  sinusoidally-oscillating  flows. 
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NOMENCLATURE 


CjL  =  in-line  force  coefficient 

Cl  =  transverse  force  coefficient 

D  =  diameter 

K  =  Keulegan-Carpenter  Number  =  U^T/D 

Ps  =  pressure  on  cylinder 

R  =  cylinder  radius 

Re  =  Reynolds  number  =  UmD/v 

r  =  radial  distance 

T  =  period  of  oscillations 

t  =  time 

U  =  time  dependent  velocity 

Um  =  maximum  velocity  in  pure  sinusoidal  flow 

p  =D2/vT 

p  =  dynamic  viscosity 

V  =  kinematic  viscosity 

p  =  density 

0  =  angular  position 

i|»  =  stream  flmction 

(i)  =  vorticity 
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1.  INTRODUCTION 


The  study  of  time-dependent  flows  past  bluff  bodies  has  historically  been  the  focus 
of  a  great  deal  of  scientific  attention  owing  to  its  relevance  to  many  and  diverse  applications, 
ranging  from  submerged  structures  to  hot  wire  anemometers.  Some  forty  years  ago, 
Keulegan  and  Carpenter  commenced  the  methodical  investigation  of  the  fluid-structure 
interactions  which  occur  when  bodies  are  immersed  in  unsteadily  flowing  fluids.  Today,  the 
effect  of  the  pertinent  parameters,  such  as  Reynolds  number,  Keulegan-Carpenter  number 
and  relative  roughness,  to  name  a  few,  is  much  better  understood  thanks  to  ongoing  research 
in  this  area.  In  spite  of  this  progress,  real  time  prediction  of  the  forces  caused  by  unsteady 
flows  on  submerged  objects,  such  as  those  acting  on  an  underwater  robotic  arm  or  an 
offshore  oil  rig,  seriously  challenges  our  current  theoretical  and  computational  capabilities. 

One  of  the  principal  empirical  tools  used  by  the  engineer  to  solve  the  problems 
described  above,  Morison's  equation,  is  only  reliable  in  predicting  forces  and  moments  in 
highly  idealized  conditions  with  very  low  or  very  high  flow  oscillation  periods. 
Additionally,  most  real-life  problems  involve  flows  in  the  turbulent  regime,  further 
increasing  the  level  of  difficulty  of  both  analytical  and  numerical  solutions. 

Given  the  cost  and  complexity  of  the  laboratory  apparatuses  required  to  reproduce 
real  life  flow  conditions,  theoretical  and  experimental  advances  have  been  paralleled  by  a 
considerable  research  effort  in  the  area  of  computational  fluid  dynamics  to  assist  in  the 
solution  of  problems  related  to  bodies  immersed  in  imsteadily  flowing  fluids.  Progress  in 
numerical  techniques  and  the  ever  increasing  power  of  present  day  computers  are  finally 
making  it  possible  to  overcome  the  barrier  historically  imposed  by  the  physical  and 
numerical  instabilities  which  have  caused  the  modeling  of  turbulent  flows  particularly 
challenging  in  the  past. 

Whereas,  imtil  recently,  numerical  experimentation  has  only  been  successful  in 
determining  flow  patterns,  Strouhal  numbers  and  force  coefficients  for  flow  regimes  within 
selected  ranges  of  relatively  low  Reynolds  numbers,  it  is  now  possible  to  obtain  useful  data 
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for  more  turbulent  flows  exhibiting  extensive  separation.  Furthermore,  commercial  software 
is  becoming  available  which  affords  researchers  much  greater  flexibility  than  the  ad  hoc 
codes  previously  generated  to  solve  very  specific  categories  of  problems,  allowing  for  the 
analysis  of  a  much  broader  gamut  of  flow  conditions. 

Unfortunately,  in  spite  of  such  great  flexibility,  modem  software  does  not  yet 
absolve  the  user  from  having  to  fine-tune  the  code  by  the  judicious  selection  of  parameters, 
numerical  techniques  and  turbulence  models.  Moreover,  even  after  the  available  numerical 
schemes  have  undergone  a  considerable  amount  of  refinement,  experimentally  obtained 
results  cannot  be  exactly  duplicated  in  all  cases. 

If  complete  and  accurate  solutions  are  not  yet  always  achievable,  however, 
approximate  solutions  can  certainly  aid  in  expanding  our  current  level  of  understanding  and 
provide  the  engineer  with  a  valuable  tool  to  improve  design  optimization.  From  the  above 
discussion,  it  can  be  inferred  that  a  benchmark  body  of  reliable  calculations,  performed  using 
the  most  sophisticated  and  carefully  selected  set  of  computational  tools,  is  required  to 
calibrate  our  current  empirical  equations  and  experimental  results. 

The  objective  of  this  investigation  is  to  lay  the  foundations  for  such  an  endeavor  and 
improve  our  current  ability  to  predict  the  effects  of  time-dependent  turbulent  flows  over 
circular  cylinders  through  the  use  of  a  state  of  the  art  commercial  code  generated  by  CFD 
Research  Corporation  (CFDRC).  Whenever  possible,  the  results  achieved  have  been 
compared  to  those  obtained  experimentally.  Besides  their  intrinsic  value,  these  results  will 
aid  in  pointing  out  some  of  the  strengths  and  weaknesses  of  the  code  and  hopefully  add  to 
our  understanding  of  the  physics  of  flow  types  not  previously  analyzed. 
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II.  BACKGROUND  STUDIES 


Research  at  the  Naval  Postgraduate  School  (Sarpkaya  1976)  resulted  in  a 
comprehensive  series  of  experiments  concerning  sinusoidally  oscillating  flow  about  smooth 
and  rough  bodies,  and  introduced  the  parameter  p  (=  Re/K)  to  assess  the  influence  of  scale 
in  periodic  flows.  The  findings  demonstrated  a  dependence  of  the  force-transfer  coefficients 
on  Re,  K,  and  k/D.  The  work  continued  involving  a  coexisting  flow  as  specified  by  U  =  Uq 
+  Un,  sin  (27tt/T)  in  which  Uo  is  the  steady  mean  velocity  and  is  the  amplitude  of 
sinusoidal  oscillations  (Storm  1980). 

Numerous  numerical  predictions  of  the  Strouhal  number,  the  pressure  distribution, 
and  the  evolution  of  the  lift  and  drag  forces  in  impulsively-started  steady  ambient  flows  have 
been  performed  to  replicate  this  effort.  No  computational  attempt  was  made  to  investigate 
the  effect  of  the  initial  acceleration,  prior  to  the  establishment  of  a  steady  uniform  flow,  on 
the  characteristics  of  the  resulting  time-dependent  flow  (Putzig  1991).  Here,  only  the  more 
recent  investigations  will  be  reviewed. 

A  finite-difference  analysis  of  the  Navier-Stokes  equations  for  a  sinusoidally 
oscillating  ambient  flow  about  a  circular  cylinder  at  K  =  5  (Re  =  1000)  and  K  =  7  (Re  =  700) 
has  been  attempted  (Baba  Miyata  1987),  assuming  a  physically  unrealistic  symmetric  wake 
in  both  simulations.  Their  results  have  shown  that  the  calculations  can  be  carried  out  only 
for  short  times  (less  than  two  cycles  of  flow  oscillation)  with  a  non-super  computer. 

A  similar  method  was  used  to  analyze  three  cases  (K  =  5,  7,  and  10)  at  higher 
Reynolds  numbers  around  1 0^  (Murashige  1 989).  The  flow  was  perturbed  by  artificial  means 
to  trigger  asymmetry.  At  K=  10,  a  transverse  vortex  street  appeared,  in  agreement  with 
experimental  observations.  A  Simulation  using  multi-discrete  vortices  (with  wore), 
simulated  the  sinusoidally  oscillating  flow  about  a  circular  cylinder  and  the  decelerating  flow 
about  camber  plates  (Mostafa  1987).  Hjs  calculations  for  K  =  12  reproduced  correctly  for 
the  first  time  the  transverse  vortex  street  observed  experimentally.  However,  the  calculated 
forces  were  somewhat  larger  than  those  measured  (Putzig  1991). 
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A  series  of  numerical  experiments  with  sinusoidally  oscillating  as  well  as  co-existing 
laminar  flows  at  relatively  small  K  values  yielded  total  force  coefficients  which  were  not 
only  highly  stable,  but  also  in  reasonable  agreement  with  those  obtained  experimentally.  The 
numerical  experiments  with  co-existing  flows  produced  extremely  interesting  flow  features. 
For  relative  current  velocities  V,  (UyU^  =  0.7  -  0.8,  the  vortices  shed  nearly  symmetrically 
at  each  cycle.  However,  they  formed  a  three-row  vortex  street  only  in  the  range  of  =  0.6- 
0.7.  For  V,  larger  than  about  one,  the  vortex  wake  returned  to  the  asymmetric  mode,  as  in 
a  regular  vortex  street.  The  calculations  of  resistance  in  co-existing  flows  showed  that  the 
inertia  as  well  as  the  drag  coefficients  for  K  =  4  -  6  were  in  reasonable  agreement  with 
experimental  data  (Sarpkaya,  Putzig,  Gordon,  Wang,  &  Dalton  1992). 

Work  has  also  centered  on  the  square  cylinder.  Early  work  was  flawed  by  the  use  of 
central  differencing  at  large  cell  Reynolds  numbers.  This  led  to  spatial  oscillations  ahead 
of  the  rectangle.  Improvements  were  made  by  using  a  time  differencing  and  third-order 
upwinding  on  the  convective  terms,  although,  because  of  standard  centered-diffusion 
differencing,  it  was  overall  second-order  accurate  spatially  for  non-zero  kinematic  viscosity. 
This  was  tried  for  a  Reynolds  number  under  3,000  (Davis  &  Moore  1982).  The  use  of  the 
modified  K-epsilon  (k-e)  and  K-omega  (k-Q)  models  on  a  square  cylinder  led  to  agreement 
comparable  with  that  obtained  formerly  only  with  the  second-moment  closure,  and  with 
respect  to  the  turbulence  energy  led  to  markedly  superior  behavior  in  the  near-field  region 
(Kato  &  Launder  1993). 

Aerodynamic  research  on  airfoils  further  showed  the  capabilities  of  the  k-e  model. 
The  model  did  require  more  grid  points  than  previous  models  like  the  Baldwin-Barth  model 
but  was  significantly  better  at  computing  maximum  lift  conditions  and  flap  boundary-layer 
separation  (Rogers  1994). 

Most  of  the  initial  numerical  codes  were  university  or  government  (NASA)  generated 
but  now  the  commercial  sector  has  started  to  produced  software  for  fluid  flows  and  data 
visualization.  The  CFD  Research  Corporation  (CFDRC)  created  a  general-purpose 
Computational  Fluid  Dynamics  (CFD)  code  with  multi-domain  solution  capability.  The 
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initial  results  showed  the  program  to  be  able  to  predict  Strouhal  numbers  in  uniform  flow  at 
higher  Reynolds  numbers  and  forces  but  questioned  the  ability  to  capture  the  high  turbulence 
intensity  levels  present  in  the  near-wake  region  (Singhal  &  Awa  1994).  The  program  was 
further  tested  in  uniform  flow  using  the  k-e  and  RNG  (Re-Normalization  Group)  model  at 
Reynolds  numbers  over  10®.  It  was  able  to  reasonably  reproduce  Strouhal  numbers  but  was 
not  tested  for  force  coefficients  (Habchi  &  Hufford  1995). 
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III.  NUMERICAL  REPRESENTATION 


A.  COMPUTATIONAL  METHOD 

Here  only  a  brief  description  of  the  computational  method  is  presented.  A  more  in 
depth  description  is  given  by  CFD-ACE  Theory  Manual. 

The  fluid  is  assumed  to  be  two-dimensional,  incompressible  and  viscous.  The 
governing  equations  of  the  code  are  the  Favre, (density)  Averaged  Navier-Stokes  (FANS) 
equations  which  are  solved  by  a  pressure- velocity  coupled  method. 

The  Navier-Stokes  equations  are  averaged  over  time  or  ensemble  of  statistically 
equivalent  flows  to  yield  averaged  equations.  In  the  averaging  process,  a  flow  quantity 
is  decomposed  into  mean  and  fluctuating  parts.  The  Favre  averages  use 


u  =  u  *  ul'  where  «=  p<j)/p 

The  overbar  denotes  Reynolds  averaging  while  tilde  denotes  Favre  averaging.  The  time 
period  of  averaging,  T,  should  be  large  compared  to  the  fluctuation  time  scale  so  that  mean 
quantities  are  stationary  over  a  number  of  samples. 

The  Favre  averaged  continuity  equation  is 


dp 

dt 


_d_ 

dx. 

J 


(pu)  -  0 


(2) 


Similarly  when  the  Navier-Stokes  (NS)  equations  are  averaged,  the  following  FANS 
equation  is  obtained  (Cebeci  &  Smith  1974): 


dt 


dx^.  '  dx. 


(3) 


The  FANS  equation  contains  less  information  than  the  full  NS  equations,  but  have  additional 
unknown  terms  -  called  the  Reynolds  stresses,  (■?“/“>).  The  correlations  between  the 
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fluctuation  components  arise  in  the  averaging  process  and  require  modeling  for  closure  of 
the  FANS  equation.  The  program  treats  the  Reynolds  stress  as  a  linear  function  of  the  mean 
strain  and  is  incorporated  in  all  the  turbulence  models  of  the  code. 

Each  flow  variable  is  governed  by  a  partial  differential  equation  (PDE)  which  is 
numerically  solved  to  obtain  a  discrete  solution  for  that  variable.  No  governing  PDE  for 
pressure  is  present.  Pressure-based  methods  utilize  the  continuity  equation  to  formulate  an 
equation  for  pressure.  The  code  uses  a  variant  of  the  Semi-Implicit  Method  for  Pressure 
Linked  Equations  (SIMPLE)  (Comini  &  Del  Giudice  1987)  which  iteratively  creates  an 
equation  for  pressure-correction  from  the  continuity  equation.  The  key  in  the  calculation  of 
the  velocity  field  is  the  unknown  pressure  gradients.  The  pressure  gradients  can  be  written 
as  sums  of  estimated  (*)  and  correction  (')  values: 

dy'^dy’  ^  (4) 


where  the  pressure  gradients  at  the  previous  step  n  make  an  estimate.  The  conservation  of 
momentum  can  be  linearized  and  split  for  a  constant  property  and  two-dimensional  flow  to 
yield  for  the  v  case: 


^dt  dt  Re  dy  ^  By 


d  ,  9v ,  „dv' 

— (P.— )]-(v"-— 

dz  dz  dy 


dz  dy  dy 


(5) 


Solving  for  the  components  of  accelerations  marked  with  asterisks  via  an  implicit  scheme 
in  the  time  integration: 


=  A(p  ^)]-(v 

dt  Re  dy  ‘  dy  dz  ‘  dz  dy  dz  dy 


(6) 


and 


c^v-Ip!  (”7) 

dt  dy 
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After  enforcing  continuity  and  rearranging,  the  following  is  obtained: 

EFp'  SFp'  1  ,5v‘ 
dy^  dz^  At  dy  dz 


(8) 


Equation  (8)  is  a  Poisson-like  equation  for  the  pressure  correction  whose  initial 
boundary  conditions  are: 

p'  =  0  (9) 

where  the  pressure  is  known  and: 


(10) 


elsewhere  .  Once  the  pressure  correction  has  been  determined,  the  corresponding  velocity 
corrections  v'  can  be  computed  from  equations  (7).  The  pressures  and  velocities  can  be 
updated: 


pn+l  =pn  +  pi 

(11) 

V"*'=v‘--^^At  =  V‘  +  v' 

(12) 

dy 

The  important  operations,  in  the  order  of  execution,  are: 

1 .  Estimate  (guess)  the  pressure  field; 

2.  Using  the  estimated  pressure,  solve  the  momentum  equations  to  obtain 
approximate  velocity  v*; 

3.  By  equation  (8)  calculate  the  pressure  correction  p'  that  enforces  continuity; 

4.  Use  the  pressure  correction  p'  to  update  the  pressure  and  velocity  fields; 

5.  Return  to  step  2)  and  repeat  until  convergence,  or  steady  state,  has  been 
reached. 
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B.  CALCULATION  OF  THE  FORCE  COEFFICIENTS 

The  in-line  and  transverse  force  coefficients  are  determined  from  the  combined 
contributions  of  the  shear  and  pressure  forces  acting  on  the  cylinder.  The  viscous  forces  are 
calculated  from  ~  M'"-  The  total  in-line  force  then  reduces  to 

2it  2n 

^iL  =  -  j(p,*cos(d)Rdd  -  Jpa)sm(0)c/0 

0  0 

and  the  total  lift  force  as 

2ir  2it 

=  -  J(p^*cos(Q)Jtdd  -  JfiG)sm(0)J0 

The  program  performs  this  integration  by  calculating  the  shear  and  pressure  forces 
on  the  wall  for  each  Cartesian  axis.  The  pressure  in  each  cell  block,  which  has  a  wall 
boimdary,  is  subdivided  into  its  coordinate  components  and  then  multiplied  by  its  area 
normal  to  that  direction.  The  pressure  forces  are  then  summed  for  all  of  the  cell  block  wall 
boundaries.  The  shear  force  is  calculated  in  each  cell  by: 


F^Shear  = 


where  6  is  the  distance  from  the  cell  center  to  the  cell  wall,  A,^  is  the  cell  wall  length 

tangential  to  a  given  direction  and  the  u  is  the  coordinate  velocity  at  the  cell  center.  The  total 

shear  force  in  a  given  direction  is  then  found  by  summing  the  cell  forces. 

The  force  coefficients  are  then  formulated: 

EF  Pressure  *11F  Shear 
C  =  _ - _ _ _ 

""  O.SpU^D  n6) 


EF  Pressure +T,F  Shear 

Q  ^  y _ y _ 

’  O.SpU^Z) 
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C.  GRID  DEFINITION 

The  choice  of  grid  definition  is  significant  in  capturing  the  gradients  of  both  the 
vorticity  and  the  stream  function,  particularly  where  they  are  large.  The  region  near  the 
cylinder  wall  is  the  most  critical  one  since  this  is  where  the  relative  velocity  gradients 
necessitates  a  higher  density  of  cell  mesh.  This  same  density  can  be  extended  radially  to 
the  desired  boundary  size  to  maintain  the  accuracy  of  the  vortices  generated  beyond  the 
separation  point.  This,  in  turn  requires  a  high  number  of  cells  which  increases  CPU  time. 
To  obtain  a  solution  in  a  realistic  time  of  24  to  48  hours  for  just  a  couple  of  vortex  shedding 
events,  one  must  make  a  trade  off  between  the  desire  to  fully  capture  the  fluid  dynamics  and 
the  computational  time.  The  most  common  way  is  to  limit  the  grid  density  in  less  critical 
regions.  The  penalty  is  the  loss  in  accuracy  of  vortex  strengths  as  they  move  away  from  the 
wall.  In  steady  flow,  this  is  not  as  critical.  A  vortex  influences  the  upstream  flow  but  does 
not  return  to  the  cylinder.  However,  this  is  not  true  in  oscillating  flow  and  the  judicial 
consideration  of  the  wake-return  problem  is  of  utmost  importance. 

A  square  grid  of  uniform  density  can  easily  be  generated  for  a  uniform  flow  about 
a  square  cylinder.  Such  a  grid  will  capture  any  event,  at  any  point,  but  computationally  it 
is  not  very  efficient  in  the  areas  away  from  the  cylinder  where  less  grid  density  is  desired. 
The  symmetry  of  the  geometry  does  reduce  the  grid  generation  and  CPU  time.  The  square 
is  not  the  primary  concern  of  this  investigation  and  a  purely  square  mesh  does  not  work  for 
circular  cylinders. 

A  polar  plot,  with  the  circular  cylinder  in  the  center,  with  radially  increasing  cell  size 
is  the  optimum  solution,  but  the  boundaries  for  flow  limits  the  viability  of  this  grid.  A  cross 
shape  quadrant  approach  can  still  capture  in  most  of  the  desired  effects  as  shown  in  Figure 
1 .  To  achieve  a  higher  density  of  mesh  points  near  the  cylinder  surface,  a  power  distribution 
may  be  utilized  in  the  radial  direction.  The  actual  specifics  of  the  grid  is  described  in 
Chapter  IV. 
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D.  BOUNDARY  CONDITIONS 


Four  boundary  conditions  are  utilized  in  this  investigation:  symmetry,  wall, 
prescribed  pressure  exit,  and  interface.  At  a  symmetry  boundary,  the  velocity  normal  to  the 
boundary  is  set  to  zero  and,  for  all  variables,  the  gradient  normal  to  the  boundary  is  set  to 
zero.  The  velocity  normal  to  the  boundary  is  also  set  to  zero  at  wall  boundaries.  In  both 
cases,  this  results  in  no  mass  source  for  the  pressure  correction  equation.  The  boundary 
pressure  correction  value  of  zero  is  linked  to  the  pressure  correction  equation  in  the 
prescribed  pressure  exit  which  results  in  mass  flow  into  or  out  of  the  region.  Interface  is 
used  for  continuity  between  domains  (regions)  of  the  grid  and  has  no  property  definition 
capabilities. 

The  program  has  two  additional  types  of  exit  and  inlet  boundary  conditions,  but  these 
are  not  used.  The  exit  boxmdary  conditions  are  an  extrapolation  intended  for  use  in 
supersonic  flow  models.  The  inlet  boundary  conditions  include  prescribed  mass  flux  and 
prescribed  total  pressure.  Inlet  boundaries,  as  defined  by  the  program,  allow  only  in  flux  of 
mass.  Thus,  to  simulate  boimdaries  that  are  subjected  to  by-directional  flow,  outlet 
prescribed  pressure  boundaries  are  utilized  (CFDRC  1993). 

E.  FLOW  GENERATION 

The  desired  flow  for  the  investigation  is  an  oscillating  flow  specified  by: 


C/  =  C/  *  C/„sin(27tt/r) 


(18) 


where  Up  is  the  steady  mean  velocity,  T  is  the  period  of  the  oscillation,  t  is  time  and,  Un,  is 
the  amplitude  of  sinusoidal  oscillations.  The  code  allows  entries  of  flow  by  velocity. 
However,  oscillating  flow  could  not  be  generated  in  this  manner.  An  alternate  manner  of 
control  by  boundary  pressure  does  generate  the  desired  velocity.  By  placing  equal 
magnitude  positive  and  negative  pressures  on  opposite  sides  of  the  boundaries,  the  required 
velocity  can  be  formulated  by: 
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(19) 


Pressure  =  —  oU  hx 

ijy  •  M 

where  ax  is  the  distance  from  the  boundary  wall  to  the  center  of  the  cylinder  for  a  symmetric 
mesh. 

D.  PARAMETER  OPTIONS 

The  CFD-ACE  is  a  menu  driven  program  that  creates  a  formatted  input  file  which  can 
be  externally  edited.  This  section  will  discuss  only  those  parameters  that  apply  to  exterior 
flow  on  bluff  bodies  and  require  additional  emphasis  beyond  the  CFD-Command  Language 
Manual. 

The  program  has  two  temporal  options.  The  default  is  Euler  and  works  with  no 
difficulties.  Crank-Nicholson  is  the  non-default  variable.  The  menu  creates  the  option  but 
fails  to  generate  the  correct  spelling  and  requires  exterior  editing  of  the  input  file  by  changing 
the  command  to  "CRANK_NICHOLSON". 

Four  spatial  options  are  available  for  the  user.  The  most  robust  is  the  default  upwind, 
which  is  satisfactory  for  most  problems.  The  second  order  and  Osher  options  fail  to  work 
on  NFS’s  SGI  R8000  but  this  is  being  corrected  with  an  update  at  the  time  of  this  writing. 
The  last  option  is  central  differencing.  The  option  diverges  when  used  with  finely  defined 
meshes  but  the  program  allows  the  user  to  modify  the  degree  of  incorporation  with  the 
upwind.  A  0.3  initial  value  of  blending  works  with  most  grid  patterns,  while  complex  three 
dimensional  problems  require  a  more  robust  0.5  mixture.  The  user  manual  does  not  give  an 
example.  The  following  example  command  for  a  turbulent  case  is  offered  to  assist  follow 
on  researchers: 

S_SCHEME  UPWIND  RHO  K  D 

S_SCHEME  CENTRAL  U  V 

S_BLENDING  0.3  U  V 

All  input  properties  are  in  metric  dimensions.  The  time  is  specified  by  a  start  and 
stop  time  in  seconds.  The  time  increment  for  steps  is  specified  by  imputing  the  number  of 
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twenty  steps  will  yield  a  0.1  second  step  increment.  An  eighty  steps  per  vortex  shedding 
period  is  a  rough  refinement  for  steady  flow  simulations. 

The  program  has  the  ability  to  input  polynomial  and  sine  or  cosine  functions,  which 
uses  radians  for  their  units  vice  degrees,  for  inputs.  The  program  limits  combinations  to 
only  two  functions  for  any  single  parameter.  The  use  of  the  spline  feature  is  a  substitute 
manner  of  input  which  allows  the  user  to  create  more  complete  situations  but  it  is  limited  to 
one  himdred  input  points.  By  using  the  restart  feature,  a  user  can  sequentially  run  any 
desired  length  of  problem  and  effectively  create  an  endless  spline.  Testing  shows  no  loss  of 
accuracy  in  computations  by  using  a  restarted  spline. 

The  code  allows  the  operator  several  options  for  output.  By  using  the  "PRINTF 
WALLS"  command  in  the  output  section,  the  Yplus  values  of  the  wall  will  be  listed  in  the 
output  file.  A  value  of  approximately  eleven  or  less  indicates  the  Low-Re  model  should  be 
used  for  turbulent  models.  Through  the  use  of  external  Fortran  programs,  a  user  can 
summarize  and  collect  the  force  data  for  additional  analysis.  This  investigation  uses  Fortran 
and  MATLAB  programs  to  produce  the  force  coefficient  curves  found  in  the  Appendix. 

The  computer  code  is  mainly  used  for  other  applications  where  the  minimum  relative 
pressure  is  not  negative.  This  is  applicable  for  interior  flows,  however,  this  is  not  true  for 
exterior  flows  on  bluff  bodies.  The  minimum  pressure  must  be  modified  in  the  output 
section  of  the  program  to  include  the  following  command  to  obtain  correct  results: 

P  FORCE  ON  PMIN  = -lE+05 
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IV.  DISCUSSION  OF  RESULTS 


A.  INTRODUCTION 

The  numerical  experiments  are  carried  out  through  the  use  of  a  Silicon  Graphics 
R8000  computer  with  a  two  gigabyte  hard  drive.  The  investigation  uses  the  CFD-GEOM 
program  as  the  grid  generator  and  the  CFD-ACE  program  as  the  equation  solver.  The  code 
is  relatively  new  and  has  very  few  papers  describing  the  results  of  previous  work.  Therefore, 
the  investigation  starts  with  low  steady  flow  (Re  <  500)  trails  to  verify  the  laminar 
capabilities.  This  is  followed  by  a  high  flow  case  (Re  =  5  x  10'*)  to  test  the  turbulent  features 
of  the  program.  An  oscillating  flow  case  of  p  =  1,000  and  K  =  20  is  investigated  for  its 
sensitivity  to  variable  code  parameters  followed  by  a  comparison  of  various  p  and  K  cases 
with  historical  laboratory  experiments. 

B.  LOW  FLOW  VERIFICATION 

Low  flow  verification  is  performed  on  two  types  of  cylinders;  square  and  round.  The 
first  testing  is  perform  on  a  square  cylinder  at  Re  =  400  for  the  purpose  of  familiarization 
and  visualization.  A  0.01  m  block  is  generated  on  a  grid  and  subjected  to  a  0.04  m/s  uniform 
steady  flow.  The  results  show  a  Karman  Street  with  no  imputed  disturbance  as  seen  in 
Figure  2. 

The  round  cylinder  is  subjected  to  the  same  conditions  but  it  generates  two  symmetric 
downstream  vortices  as  seen  in  Figure  3.  The  vortices  become  1  diameter  in  length  at 
approximately  1 .5  diameters  of  flow  and  then  expand  lengthwise  and  in  magnitude  for  the 
remainder  of  the  testing.  This  result  matches  laboratory  experiments  for  non-disturbed  cases. 
If  perturbed,  a  Karman  Street  forms.  The  geometry  of  the  square  cylinder  creates  enough 
discontinuities  at  the  comers  so  that  no  artificial  disturbance  is  necessary  for  that  case. 

A  simulated  disturbance  of  a  ramped  25  %  increase  and  decrease  in  steady  flow  rate 
at  the  top  half  of  the  cylinder  and  proportionally  the  same  amount  of  decrease  and  increase 
at  the  bottom  half  is  added  to  create  a  cross  flow  on  the  formed  vortices.  It  is  timed  to  start 
at  1 .5  diameters  of  flow  for  1  diameter  of  flow.  The  lower  the  percentage,  the  longer  the 
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disturbance  time  has  to  run  to  obtain  the  desired  cross  flow  condition.  This  creates  a 
Karman  Street,  as  shown  in  Figure  4.  Table  1  shows  a  comparison  of  the  experimental  and 
historical  results.  The  Strouhal  number  and  the  force  coefficient  do  not  compare  favorably 
at  Re  =  100.  However,  the  testing  does  show  an  improving  trend  for  both,  with  increasing 
speeds  of  flow.  This  test  indicates  that  the  model  does  not  simulate  properly  the  initiation 
and  magnitude  of  vortices  in  the  near  wall  region,  but  as  the  flow  increases,  which  lessens 
the  importance  of  these  effects,  the  results  improve. 

Table  1.  Comparison  of  Experimental  and  Historical  In-Line  Force  Coefficients 

and  Strouhal  Numbers  at  Low  Reynolds  Numbers 


CjL  Strouhal  Number 


Re 

Experimental 

Historical 

Experimental 

Historical 

100 

1.07 

1.8 

0.13 

0.165 

200 

1.0 

1.4 

0.14 

0.18 

400 

1.1 

1.2 

0.16 

0.19 

C.  HIGH  FLOW  VERIFICATION 

Previous  testing  (Hufford  1995)  shows  the  model  to  be  reasonably  accurate  for 
steady  flow  about  a  square  cylinder.  However,  the  area  of  interest  for  this  investigation  is 
a  circular  cylinder.  A  single  test  at  Re  =  5x10“  is  conducted  with  a  turbulence  level  of  3 
percent,  Euler  temporal  differencing,  central  spacial  differencing  and,  a  1/400  time  interval 
to  try  to  obtain  approximately  80  steps  per  vortex  separation.  This  yielded  a  Cjl  =  0.95 
which  compares  favorably  to  the  historical  value  of  1 .2.  For  two  cycles,  a  Strouhal  number 
of  0.24  is  obtained  which  compares  to  a  historical  value  of  0.23.  Based  on  this  result  and  the 
low  Re  testing,  the  model  is  found  to  sufficiently  accurate  to  examine  oscillatory  flows  at 
higher  Re  numbers. 
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D.  OSCILLATING  FLOWS 

The  baseline  grid  used  for  all  of  the  testing  of  oscillating  flow  is  a  1  meter  by  1  meter 
square  with  a  0.05  meter  diameter  circular  cylinder  positioned  in  the  center.  There  are  65 
equally  spaced  mesh  points  on  each  of  the  outer  boundaries  while  130  points,  with  a  .975 
exponential  spacing,  are  located  in  the  axially  direction.  This  creates  8,256  cells  in  each 
quadrant  for  a  total  of  33,024  for  the  model.  This  combination  is  used  to  generate  cells  that 
are  half  as  long  in  the  radial  direction  as  their  length  along  the  cylinder  wall  while  achieving 
nearly  square  cells  near  the  boundary  to  better  capture  the  pressure  and  velocity  gradients  by 
the  wall.  The  resulting  structure  can  be  seen  in  Figure  1 . 

The  initial  trial  run  uses  a  K  -  e  turbulence  model  and  SIMPLEC  algorithm.  The  time 
interval  is  0.02  seconds  with  a  first  order  upwind  spatial  differencing  scheme  and  euler 
temporal  differencing  scheme.  The  turbulence  level  at  the  boundary  is  3  %  with  a 
corresponding  turbulence  length  scale.  No  disturbances  are  required  to  generate  the  lift  crisis 
for  oscillating  flow  cases  and  therefore  not  used.  Twenty  iterations  are  performed  at  each 
time  step.  All  runs  use  p=  1,000  and  K  =  20  for  a  Re  =  20,000.  This  flow  velocity 
combined  with  the  grid  geometry  and  time  intervals  creates  a  Neumann  number  of  0.5  at  the 
boundaries  of  the  grid  area. 

The  code's  variable  parameters  are  tested  for  their  sensitivity  and  the  results  are 
summarized  in  Table  2.  The  initial  baseline  model  obtains  a  Cq,  of  1.1  compared  to  a 
historical  maximum  peak  value  of  1.9. 
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Table  2.  Sensitivity  of  Parameters  on  the  Cn, 


Relative  Percent 


Parameter 

Factor  of  Change 

of  Change 

Original 

New 

Time  Interval 

0.02  to  .01 

0.5 

1.33 

1.34 

Time  Interval 

0.02  to  .0025 

0.5 

1.33 

1.34 

Turbulence  Level 

3%  to  10% 

2.4 

1.80 

1.85 

Model  Type 

K  -  e  to  RNG 

6.8 

1.09 

1.22 

Model  Type 

K  -  e  to  Low  Re 

12.3 

1.09 

1.33 

Spatial  Differencing 

Upwind  to  Central 

9.7 

1.33 

1.52 

Temporal  Differencing 

Euler  to 

-4.1 

1.33 

1.25 

Crank  Nicholson 

Turbulence  Length  Scale  0.0025  to  0.0045 

2.0 

1.09 

1.13 

Number  of  Iterations 

40  to  20 

0.0 

1.22 

1.22 

per  step 

Number  of  Iterations  Step  20  to  10 

0.0 

1.22 

1.22 

per  Step 

Number  of  Iterations 

20  to  5 

-5.1 

1.22 

1.22 

per  Step 

The  optimum  combination  uses  10  iterations  per  step  with  a  0.02  time  interval  per 
step  which  does  not  impact  the  results  but  reduces  the  time  to  run  the  problem  by  over  60  % 
compared  to  the  original  model.  The  temporal  differencing,  turbulence  level  and,  length 
scale  all  remain  the  same  as  the  baseline  model.  The  simulation  is  altered  to  have  a  Low  Re 
turbulence  model  and  central  spatial  differencing.  This  combination  increases  the 
maximum  observed  Cn,  to  1.5.  Additional  testing  shows  the  lift  crisis  significantly  increases 
the  peak  Cq,  from  1 .5  to  1 .8  which  may  of  contributed  to  the  low  value  for  the  initial  model. 
Figure  5  shows  the  visualization  of  one  cycle  of  flow. 
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The  grid  size  definition  is  tested  next  by  doubling  the  number  of  boundary  points 
creating  a  133,644  cell  grid.  The  grid  refinement  is  able  to  capture  more  of  the  magnitude 
of  the  vortices  than  the  course  grid  which  results  in  a  1 .7  Cn,  prior  to  lift  crisis  and  a  1 .8  after 
the  crises.  However,  this  increase  in  grid  definition  slows  the  computer  by  over  400  % 
compared  to  the  basic  line  model  and  does  not  provide  additional  improvement  in  force 
coefficients  after  the  lift  crisis. 

The  optimum  model  is  next  tested  at  various  p  and  K  combinations  that  correspond 
to  previous  testing  for  comparison  (Sarpkaya  1976).  Figure  5-34  show  the  in-line  and 
transverse  forces  for  each  combination  and  are  summarized  in  Table  3. 

Table  3.  Oscillating  Flow  In-Line  Force  Coefficient  Results 


K 

p  =  1,000 

Computer  Historical 

p  =  2,000 

Computer  Historical 

p  =  3,100 

Computer  Historical 

8 

2.6 

2.4 

2.45 

2.5 

2.45 

2.4 

12 

2.25 

2.3 

1.72 

2 

1.7 

1.45 

20 

1.8 

1.95 

1.15 

1.3 

1.1 

1.0 

25 

1.3 

1.6 

0.95 

1.15 

1.0 

0.9 

35 

1 

1.35 

0.83 

0.85 

0.9 

0.85 

The  following  observations  are  made; 

1 .  The  p  =  2,000  data  acted  more  like  a  higher  p  number  and  particularly  in  the 
mid  K  flows. 

2.  The  K=8  and  K=35  cases  usually  yielded  best  similarity  to  the  historical  data. 
The  K=8  case  may  be  due  to  the  flow  being  less  complex  than  K=12  case. 
At  high  Ks,  a  plateau  is  reached  where  the  Cjl  is  more  relatively  constant. 
The  code  is  able  to  capture  most  of  the  vortex  strength  and,  therefore,  the 
results  tend  to  be  more  accurate  in  the  region  of  higher  K  values. 
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3.  The  dissipation  of  visualized  vortices  and  their  strengths  are  higher  than 
nature  or  previously  used  codes.  This  may  be  the  cause  for  the  louver  force 
coefficients. 

4.  The  p  =  1 000  has  a  longer  time  to  the  lift  crises  which  may  be  due  to  lower 
Um  values. 

A  view  of  the  vortices  during  one  oscillation  for  p  =  1,000  and  K  =  8  shows  in 
Figures  35  and  36  the  formation,  detachment,  and  the  reentry  into  the  wall  region  on  flow 
reversal.  Figures  37  through  40  show  the  start  of  a  transverse  street  for  the  p  =  1,000  and 
K  =  12  case.  Figures  41  through  44  show  the  return  of  the  normal  oscillation  pattern  for  the 
p  =  1 ,000  and  K  =  20.  Figures  45  and  46  for  p  =  2,000  and  K  =  20  begin  to  show  the  effect 
of  higher  flows  where  each  half  cycle  acts  like  an  impulsive  start  and  forms  symmetric 
vortices.  This  effect  can  be  seen  more  clearly  in  a  series  of  views  for  p  =  3,100  and  K  = 
35  shown  in  Figures  47  through  50  but  the  vortices  are  more  elongated  and  do  not  linger 
near  the  cylinder  because  of  the  higher  11^. 

Overall,  the  trend  of  the  data  is  well  predicted  for  the  pure  oscillatory  case  with  no 
steady  mean  flow.  The  program  has  limitations  in  fully  capturing  the  generated  forces  and 
dissipation  but  does  allow  a  look  at  flow  levels  not  previously  tested  on  a  non-super 
computer. 
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V.  CONCLUSIONS 


The  investigation  reported  here  warranted  the  following  conclusions: 

1.  Finite  differencing  formulations  of  the  Favre  Averaged  Navier  Stokes  equations 
(FANS)  can  reasonably  solve  a  wide  range  of  Reynolds  number  flows,  using  a 
modified  SIMPLE  method. 

2.  The  laminar  flow  model  predicts  the  trend  of  Strouhal  numbers  but  does  not 
favorably  predict  low  Reynolds  number  force  coefficients.  The  comparability  with 
historical  values  does  improve  with  higher  Reynolds  number  flows.  The  force 
coefficients  correspond  to  higher  Reynolds  number  conditions,  but  the  Strouhal 
numbers  correspond  to  the  actual  or  lower  flow  conditions. 

3.  The  numerical  experiments  with  oscillating  flows  for  p  of  1,000,  2,000  and,  3,100 
yield  In-Line  Force  Coefficients  in  good  agreement  with  those  obtained 
experimentally. 

4.  The  program  does  not  fully  capture  vortex  strengths  and  prematurely  dissipates  the 
vortices  relative  to  previous  investigations. 

5.  The  numerical  experiments  presented  herein  can  be  performed  on  a  non-super 
computer.  The  fact  that  numerical  instabilities  versus  fluid  dynamical  instabilities 
have  different  and  at  times  competing  mathematical  and  computational  demands 
necessitated  a  compromise  in  grid  definition  and  time  interval  spacing.  In  spite  of 
this,  eight  hours  of  computer  run  time  were  required  to  generate  twenty  seconds  of 
true  flow  for  most  of  the  oseillatory  cases. 
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APPENDIX 


Figure  1.  Grid  in  the  Physical  Domain 


Figure  2.  Flow  About  a  Square  Cylinder  at  Re  =  400 
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Figure  5.  Inline  Force  Coefficient  vs  t/T,  P  =  1,000,  K  =8,  Re  =  8,000 


Figure  6.  Transverse  Force  Coefficient  vs  t/T,  p=  1,000,  K=8,  Re  =  8,000 
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Transverse  Force  coefficient  cfq’  In  Line  Force  Coefficient 
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Figure  9.  Inline  Force  Coefficient  vs  t/T,  P  =  1,000,  K  =20,  Re  =  20,000 
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vr 

Figure  10.  Transverse  Force  Coefficient  vs  t/T,  P=  1,000,  K=20,  Re  =20,000 
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Transverse  Force  coefficient 


Figure  16.  Transverse  Force  Coefficie 


Figure  17.  Inline  Force  Coefficient  vs  t/T,  P  =  2,000,  K  =  12,  Re  =  24,000 


Figure  18.  Transverse  Force  Coefficient  vs  t/T,  P=  2,000,  K=  12,  Re  =  24,000 
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In  Line  Force  Coefficient 


Figure  21 .  Inline  Force  Coefficient  vs  t/T,  P  =  2,000,  K  =  25,  Re  =  50,000 


Figure  22.  Transverse  Force  Coefficient  vs  t/T,  p=  2,000,  K=  25,  Re  =  50,000 
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T ransverse  Force  coefficient 


Figure  25.  Inline  Force  Coefficient  vs  t/T,  p  =  3,100,  K  =  8,  Re  =  24,800 


Figure  26.  Transverse  Force  Coefficient  vs  t/T,  P  =  3,100,  K  =  8,  Re  =  24,800 
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)t  for  an  Oscillatory  Flow  About  a  Circular  Cy! 
(Steps  1  -  6  in  Increments  of  45  Degrees,  Left  1 


Vorticity  Plot  for  an  Oscillatory  Flow  About  a  Circular  Cylinder  at  p  =  1000,  K  =  8,  F 
(Steps  7  - 12  in  Increments  of  45  Degrees,  Left  to  Right,  Top  to  Bottom) 


Figure  37.  Vorticity  Plot  for  an  Oscillatory  Flow  About  a  Circular  Cylinder  at  p  =  1000,  K  =  12,  Re  =  12,000 

(Steps  1  -  6  in  Increments  of  22.5  Degrees,  Left  to  Right,  Top  to  Bottom) 


Vorticity  Plot  for  an  Oscillatory  Flow  About  a  Circular  Cylinder  at  p  =  1000,  K  =  12,  R( 
(Steps  7  -  12  in  Increments  of  22.5  Degrees,  Left  to  Right,  Top  to  Bottom) 


Figure  39.  Vorticity  Plot  for  an  Oscillatory  Flow  About  a  Circular  Cylinder  at  p  =  1000,  K  =  12,  Re  =  12,000 

(Steps  13  -  18  in  Increments  of  22.5  Degrees,  Left  to  Right,  Top  to  Bottom) 


Vorticity  Plot  for  an  Oscillatory  Flow  About  a  Circular  Cylinder  at  p  =  1000,  K  =  12,  Re  =  12,000 
(Steps  19  -  24  in  Increments  of  22.5  Degrees,  Left  to  Right,  Top  to  Bottom) 


Figure  41 .  Vorticity  Plot  for  an  Oscillatory  Flow  About  a  Circular  Cylinder  at  p  =  1000,  K  =  20,  Re  -  20,000 

(Steps  1  -  6  in  Increments  of  45  Degrees,  Left  to  Right,  Top  to  Bottom) 
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Figure  44.  Streamline  Plot  for  an  Oscillatory  Flow  About  a  Circular  Cylinder  at  p  =  1000,  K  =  20,  Re  =  20,000 

(Steps  7  -  12  in  Increments  of  45  Degrees,  Left  to  Right,  Top  to  Bottom) 
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Figure  47.  Vorticity  Plot  for  an  Oscillatory  Flow  About  a  Circular  Cylinder  at  p  =  3100,  K  =  35,  Re  =  108,500 

(Steps  1  -  6  in  Increments  of  45  Degrees,  Left  to  Right,  Top  to  Bottom) 
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(Steps  7  - 12  in  Increments  of  45  Degrees,  Left  to  Right,  Top  to  Bottom) 


Figure  49.  Streamline  Plot  for  an  Oscillatory  Flow  About  a  Circular  Cylinder  at  p  =  3 

(Steps  1  -  6  in  Increments  of  45  Degrees,  Left  to  Right,  Top 
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